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We  present  a  generalized  Kohn-Sham  (KS)  density  functional  theory  (DFT)  based  effective  fragment 
potential  (EFP2-DFT)  method  for  the  treatment  of  solvent  effects.  Similar  to  the  original  Hartree- 
Fock  (HF)  based  potential  with  htted  parameters  for  water  (EFPl)  and  the  generalized  HE  based 
potential  (EFP2-HF),  EFP2-DFT  includes  electrostatic,  exchange-repulsion,  polarization,  and  dis¬ 
persion  potentials,  which  are  generated  for  a  chosen  DET  functional  for  a  given  isolated  molecule. 
The  method  does  not  have  htted  parameters,  except  for  implicit  parameters  within  a  chosen  func¬ 
tional  and  the  dispersion  correction  to  the  potential.  The  electrostatic  potential  is  modeled  with  a 
multipolar  expansion  at  each  atomic  center  and  bond  midpoint  using  Stone’s  distributed  multipolar 
analysis.  The  exchange-repulsion  potential  between  two  fragments  is  composed  of  the  overlap  and 
kinetic  energy  integrals  and  the  nondiagonal  KS  matrices  in  the  localized  molecular  orbital  basis. 
The  polarization  potential  is  derived  from  the  static  molecular  polarizability.  The  dispersion  poten¬ 
tial  includes  the  intermolecular  D3  dispersion  correction  of  Grimme  et  al.  [J.  Chem.  Phys.  132, 
154104  (2010)].  The  potential  generated  from  the  CAMB3LYP  functional  has  mean  unsigned  errors 
(MUEs)  with  respect  to  results  from  coupled  cluster  singles,  doubles,  and  perturbative  triples  with 
a  complete  basis  set  limit  (CCSD(T)/CBS)  extrapolation,  of  1.7,  2.2,  2.0,  and  0.5  kcal/mol,  for  the 
S22,  water-benzene  clusters,  water  clusters,  and  /r-alkane  dimers  benchmark  sets,  respectively.  The 
corresponding  EEP2-HE  errors  for  the  respective  benchmarks  are  2.41,  3.1,  1.8,  and  2.5  kcal/mol. 
Thus,  the  new  EEP2-DET-D3  method  with  the  CAMB3LYP  functional  provides  comparable  or  im¬ 
proved  results  at  lower  computational  cost  and,  therefore,  extends  the  range  of  applicability  of  EEP2 
to  larger  system  sizes.  ©  2014  AIP  Publishing  LLC.  [http://dx.doi.Org/10.1063/l.4883488] 


I.  INTRODUCTION 

Efficient  and  accurate  methods  for  describing  intermolec¬ 
ular  interactions  are  of  considerable  interest  in  many  helds 
of  chemistry,  physics,  material,  and  biological  sciences  be¬ 
cause  accurate  hrst  principles  methods  are  computationally 
demanding  with  increasing  system  size  for  materials  in  solu¬ 
tion,  such  as  biological  systems,  colloidal  quantum  dots,  and 
other  materials  in  the  condensed  phase.  To  reduce  the  compu¬ 
tational  cost,  the  effective  fragment  potential  (EEP)  method 
was  developed  to  explicitly  model  intermolecular  interactions 
for  molecular  complexes  (e.g.,  (H20)n)  by  using  quantum  me¬ 
chanical  (QM)  properties  of  a  single  molecule.  The  estab¬ 
lished  QM  parameters  can  be  stored  for  future  use  or  repro¬ 
duced  on  the  fly  as  needed.  Thus,  EEP  potentials  are  similar 
in  spirit  to  molecular  mechanics  force  helds  that  are  gener¬ 
ated  by  using  ab  initio  data  such  as  the  SIBEA  (sum  of  inter¬ 
actions  between  fragments  ab  initio  computed).^  In  the  EEP 
methods,  molecules  within  a  complex  are  replaced  by  frozen 
molecular  fragments  with  hxed  internal  coordinates.  The  EEP 
intermolecular  interactions  are  represented  by  Coulomb,  in¬ 
duction,  and  dispersion  terms  that  are  generated  by  auxiliary 
hrst  principles  calculations.  The  present  work  focuses  on  the 
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development  of  a  second  generation  EEP  (EPP2)  based  on 
density  functional  theory  (DPT).  A  number  of  reviews  have 
been  published,  including  most  recent  EEP  developments"^  on 
the  interface  with  the  ground  and  excited  state  ah  initio  wave- 
functions. 

The  hrst  generation  of  EPPs^’^  based  on  Hartree-Pock 
(HP)  and  density  functional  theory  (DPT)  were  originally  for¬ 
mulated  and  implemented  to  model  water  (called  EPPl).^’^’^ 
The  HP  based  EPP2  extension  with  dispersion  between 
fragments^  or  HP-EPP2^  systems  is  a  more  general  and 
parameter-free  method  for  any  molecular  species. How¬ 
ever,  the  EPPl-DPT^  potential  was  htted  using  the  B3LYP 
functional,  which  lacks  long-range  and  dispersion  interac¬ 
tions.  The  EPPl-HP^  method  was  htted  to  reproduce  results 
for  HP,  which  does  not  include  the  important  electron  correla¬ 
tion  effects.  The  DPT  based  EPP2  with  dispersion  interactions 
offers  advantages  of  DPT  over  HP  at  comparable  computa¬ 
tional  cost.  However,  as  the  cost  for  evaluating  the  dynamic 
polarizability  becomes  prohibitive  for  obtaining  the  disper¬ 
sion  potential  for  large  EPP2-HP  systems,  we  are  motivated  to 
include  the  simple  D3  dispersion  potential,  which  has  been  re¬ 
cently  rehned  to  high  accuracy  and  broad  applicability.^^  The 
geometry  based  D3  dispersion  potential  is  also  compatible 
with  the  DPT-EPP2  dispersion  interaction.  The  dispersion  in¬ 
teractions  for  QM  systems  with  PEP  for  HP-EPP2  practically 
require  a  separated  treatment  for  the  QM  part  calculated  with 
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dipole  integrals  and  orbital  energies.^  In  this  work,  we  imple¬ 
mented  and  evaluated  the  accuracy  of  the  DFT-EFP2  poten¬ 
tial,  using  the  mean  unsigned  errors  (MUE)  and  mean  signed 
errors  (MSE),  relative  to  a  number  of  benchmarks  from  the 
coupled  cluster  (CC)  method  with  singles,  doubles  (CCSD) 
and  noniterative  triples  with  the  extrapolation  to  the  complete 

basis  set  limit  (CCSD(T)/CBS).i4-i8 

The  paper  is  organized  as  follows.  After  providing  a  de¬ 
scription  of  the  method  (Sec.  II)  and  computational  details 
(Sec.  Ill),  we  present  results  and  discussion  for  the  bench¬ 
marking  (Sec.  IV).  Section  V  provides  a  summary  of  our  re¬ 
sults  and  conclusions  for  the  paper. 


II.  THEORY 

Considering  a  complex  with  two  or  more  non-bonded 
interacting  molecules  (fragments),  the  quantity  of  interest  is 
its  total  interaction  energy.  The  EEP  methods  for  obtain¬ 
ing  the  interaction  energy  have  been  described  in  previous 
papers. Therefore,  we  only  briefly  review  the  theory 
in  order  to  document  our  implementation  of  the  EPP2-DPT 
method  with  the  dispersion  potential.  The  EPP2  poten¬ 
tial  with  electrostatic  Coulomb  polarization 

exchange-repulsion  and  dispersion  (V^^^P)  contribu¬ 

tions  from  an  5'  fragment  can  be  written  as 


ments  as^ 

r)  =  ^  ixl(s)fc(rrk) 

a 

^  X,y,z 

^  x,y,z 

where midfa^y  are  the  QM  electric  field,  field  gradient, 
and  field  Hessian,  respectively. 

The  multipolar  expansion  at  each  atomic  center  and  bond 
midpoint  are  done  with  Stone’s  distributed  multipolar  anal¬ 
ysis  using  Kohn-Sham  (KS)  density.  As  distributed  multi¬ 
polar  expansion  of  a  molecular  electrostatic  potential  is  only 
well  represented  at  large  intermolecular  distances,  the  point 
charge  model  of  Eq.  (2)  must  account  for  overlapping  charge 
densities  by  a  distance-dependent  cutoff  function.^’ The 
nuclei  of  a  molecule  are  not  shielded  by  their  own  electron 
density,  leading  to  a  greater  attraction  for  electron  density  as¬ 
sociated  with  other  molecules.  This  effect  is  estimated  by  an 
exponential  damping  function  fit  to  the  quantum  mechanical 
electrostatic  potential  over  a  grid  points. 


K  L 

V^wiis.  r)  =  y]  yf  “(s,  r)  +  J2  r) 

k^l  /=1 

M 


where  each  fragment  indexed  by  5’  could  have  a  number  {K, 
L,  M)  of  expansion  points;  r  is  the  electronic  coordinate. 
The  term  describes  classical  electrostatic  interactions 
between  the  charge  distributions  of  the  effective  fragments. 
The  electrostatic  energy  can  be  effectively  represented  by  dis¬ 
tributed  multipoles  at  large  intermolecular  distance  while  the 
charge  overlap  of  the  fragment’s  electron  clouds  at  the  short- 
range  is  also  taken  into  account  as  given  in  Subsection  II  A. 
The  non-classical  exchange-repulsion  energy  arising 
from  electron  exchange,  which  requires  a  QM  treatment  (Sub¬ 
section  II  B).  The  polarization  term  arises  from  the  mu¬ 
tual  perturbation  of  electrostatic  fields  of  the  unperturbed 
fragments.  The  EEP  polarization  potential  is  approximated  by 
using  polarizability  tensors  as  described  in  Subsection  II  C. 
The  dispersion  term  originates  from  induced  dipole- 

dipole  interactions  at  the  lowest  order,  which  is  treated  using 
dynamic  polarizabilities.  The  treatment  of  dispersion  for  EEP 
is  described  in  Subsection  II  D. 


A.  Electrostatic  potential 

The  electrostatic  potential  consists  of  a  multicenter  dis¬ 
tributed  multipolar  expansion  (DMA),  carried  out  with  the 
charge  {q),  dipole  (/x),  quadrupole  (©),  and  octopole  (^2)  mo- 


B.  Exchange-repulsion  potential 

The  exchange  energy  arises  from  the  Pauli  exclusion 
principle  for  electrons  with  like  spin,  which  requires  the 
wavefunctions  to  be  antisymetric  with  respect  to  electron 
exchange.  The  exchange-repulsion  potential  is  obtained  by 
subtracting  the  classical  Coulomb  energy  from  the  first-order 
intermolecular  interaction  energy  from  symmetry  adapted 
perturbation  theory,^^  which  is  identical  to  the  Heitler-London 
energy.^  By  neglecting  the  internal  energy  and  truncating  the 
intermolecular  overlap  expansion  after  second  order,  the  ap¬ 
proximated  exchange-repulsion  potential  between  two  EEP 
fragments  (A  and  B)  can  be  separated  into  three  terms: 

yEx  ^  -h  (3) 


The  zeroth  order  term  is  approximated  as 


yE’‘(5°)  =  -2V  V  (4) 

where  Kj  is  the  electron  exchange  operator.  Using  the  spheri¬ 
cal  Gaussian  overlap  (SGO)  approximation,^^  the  zeroth  order 
intermolecular  exchange  energy  for  RHP  becomes 


yEx(^O)  =  -2  y]  y]  2- 

ieA  jeB 


-2\nSij 

7t  Rij 


(5) 


where  Sij  is  the  overlap  integral  of  the  localized  molecular  or¬ 
bitals  (LMOs)  and  Rij  is  the  distance  between  LMO  centroids. 
The  first  order  term  is  approximated  as 


yHx(5i)_2y^y^5.. 

i  gA  j  gB 


Vij,A  +  Vij^B 


+  y]  {i  \2J,  -  Kk\j)  +  y]  {i  \2Ji  -  X,|  y)],  (6) 

A:gA  /gB 
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where  is  the  Coulomb  operator,  Kj^  is  the  exchange  oper¬ 
ator,  and  Vij^p,  is  the  potential  due  to  nuclei  in  A.  For  a  RHF 
wavefunction  the  first  order  term  can  be  written  in  terms  of 
the  Fock  matrixes,  overlap  integrals  (S),  and  kinetic  energy 
integrals  (T)  of  the  LMOs  as 

"A  B 

ieA  jeB  \_  k  I 

(7) 

For  DFT  based  EFP2,  F  is  the  nondiagonal  KS  matrix  in  the 
LMO  basis,  instead  of  the  Fock  matrix.  The  second  order  term 
for  the  exchange-repulsion  potential  is  approximated  as 


C.  Polarization  potential 

The  polarization  potential  is  calculated  as 


x,y,z 

=  J2  A  fe)  (s)  iffi  to)) ,  (10) 

where /is  the  field  generated  by  QM  system  and  fragments. 
The  polarizability  for  a  LMO  (j)f  by  can  be  obtained  by 
the  finite-field  method  as^^ 


a 


rs 


Fr’ 


(11) 


Skj  + 

ieA  jeB  LkeA  \  leA  / 

+  ^kj  I  (j  \^Jk  \  1)  j 

leB  \  keA  / 

+  E‘^«  + 

leB  \  keA  / 


-EE 

KeA  leB 


(8) 


The  integrals  in  (5^)  are  modeled  with  classical  electro¬ 
statics  for  both  RHF  and  DFT  as 


ieA  jeB 


E 

_7gB 


Zj  ^ 

/€B 


1 

Ft 


+E 

/eA 


keA 


(9) 


The  indexes  i,  k  and  j,  I  run  over  the  LMOs  of  fragment 
A  and  B,  respectively.  Ru  is  the  distance  between  LMO 
centroid  i  and  the  nucleus  J.  The  choice  of  DFT  function¬ 
als  affects  the  accuracy  of  EFP  exchange  repulsion  potential 
(Eqs.  (3)-(9))  that  includes  the  HE  exchange  (for  hybrid  func¬ 
tionals)  and  KS  matrices.  Although  DET  functionals  using  the 
generalized  gradient  approximation  (GGA)  provide  reason¬ 
able  energetics  for  some  weakly  bound  systems,  GGA  poten¬ 
tials  exhibit  an  exponential  decay  in  the  asymptotic  regions  of 
molecular  systems,  instead  of  the  correct  —  1/r  decay,  leading 
to  large  errors  for  intermolecular  binding  energies. Hy¬ 
brid  functionals  that  include  a  small  fraction  of  the  exact  HE 
exchange  deliver  more  accurate  results  for  thermochemistry. 
However,  these  potentials  decay  as  —dr  (where  c  is  fraction 
or  scaling  factor  of  HE  exchange),  not  the  correct  —1/r  de¬ 
cay.  Long-range  corrected  (LC)  or  Coulomb  attenuated  (CA) 
hybrids  significantly  reduce  of  the  self-interaction  errors  as¬ 
sociated  with  hybrid  functionals  by  retaining  the  full  HE  ex¬ 
change  only  for  long-range  electron-electron  interactions  in 
the  asymptotic  regions. Thus,  hybrid,  LC,  or  CA  function¬ 
als  represent  an  appropriate  choice.  The  CAMBSLYP^^  func¬ 
tional  with  available  dispersion  parameters'^  is,  therefore, 
chosen  for  the  present  work.  To  ensure  correct  long-range  in¬ 
teractions,  the  SGO  integrals  (Eq.  (5))  are  not  scaled  in  DET 
based  EEP2. 


where  /x^  is  the  5'  component  of  the  dipole  moment  induced 
by  an  electric  field  applied  along  direction  r, 

At'  =-2[(0;K|0;)-(0O|r,|0O)].  (12) 

The  induced  dipole  (Eq.  (12))  is  assumed  to  be  a  linear  func¬ 
tion  of  the  applied  field.  However,  for  the  nonuniform  fields 
encountered  in  molecular  interactions  for  systems  with  more 
than  one  fragment,  the  induced  moments  within  each  frag¬ 
ment  depend  on  the  electric  fields  of  the  other  fragment  as 
well  as  the  active  part  of  the  QM  part  of  the  system,  and 
the  fragment-fragment  polarization  energy  must  be  iterated 
to  self-consistency.  An  induced  dipole  is  placed  at  the  cen¬ 
troid  Vs  of  each  LMO  in  the  valence  shell.  Eor  example  in 
methane,  the  LMO  centroid  locations  occur  along  the  four 
C-H  bonds  axes.  Note  that  different  LMOs  result  in  different 
sets  of  polarizability  tensors.  LMO  centroid  positions  that  are 
not  well  determined  in  certain  cases  of  high  symmetry  and 
for  some  aromatic  systems,  lead  to  unphysically  large  com¬ 
ponents  of  the  polarizability  tensors.  To  obtain  the  complete 
polarizability  matrix,  at  least  three  applied  Cartesian  {x,  y,  z) 
fields  SCE  calculations  must  be  carried  out,  in  addition  to  the 
zero-field  SCE  calculations.  The  SCE  convergence  of  10“^^ 
(about  2  orders  of  magnitude  tighter  than  is  acceptable  for 
other  applications)  leads  to  induced  dipole  moments  within 
10“^  to  10“^  au  of  the  precise  values.  Thus,  the  applied  field 
of  10“"^  au  leads  to  an  uncertainty  of  10“^  to  10“^  for  the  po¬ 
larizability  using  the  equations  above.  In  the  cases  in  which 
the  finite-field  method  described  above  resulted  in  unphysi¬ 
cally  large  components  of  the  polarizability,  analytical  meth¬ 
ods  have  been  developed  using  coupled  perturbed  Hartree- 
Eock  (CPHE)  equations^^  with  exact^^“^^  and  approximate^"^ 
(A-CPHE)  response  vectors.  The  approximate  response  vec¬ 
tors  include  only  HP  orbital  energies  and  dipole  moments. 
Thus,  the  coupled  perturbed  Kohn-Sham  (CPKS)^^’^^  equa¬ 
tions  with  approximate  (A-CPKS)  response  vectors  that  in¬ 
clude  KS  orbital  energies  and  dipole  moments,  are  in  analogy 
to  the  A-CPHP  equations.^"^ 

In  the  EPP2  model,  the  polarization  potential  (Eq.  (10)) 
currently  includes  only  dipole  polarizabilities.  The  method, 
in  principle,  can  be  extended  to  any  order  of  polarization 
but  with  at  the  added  cost  and  complexity  of  analytical  gra¬ 
dients  of  polarization  energy  with  respect  to  translation  and 
rotation  motions  of  the  fragments.  However,  benchmarks  of 
symmetry-adapted  perturbation  theory  calculations  indicated 
that  the  main  source  of  EPP2  errors  come  from  the  under¬ 
estimation  of  the  Coulomb  and  polarization  terms  that  are 
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partially  compensated  by  the  underestimation  of  the  exchange 
repulsion  term.^^ 

D.  Dispersion  potential 

Calculations  of  the  dispersion  potential  for  EFP2-HF 
were  done  using  anisotropic  dynamic  polarizability  tensors 
evaluated  over  imaginary  frequencies  for  the  term  while 
higher  order  terms  were  estimated  as  1/3  of  the  term.  The 
method  has  been  described  in  detail  previously.^’ ^  The  D3  dis¬ 
persion  potential  is  calculated  as^^’^^ 

^AB 

=  E  E  (13) 

AB  n=6,8,.. 

where  C^^are  averaged  isotropic  dispersion  coefficients  for 
AB  atomic  pairs,  xab  is  the  internuclear  distance.  Thus,  in¬ 
tramolecular  dispersion  interactions  are  not  included  for  the 
fragments.  The  scaling  factors  Sn  are  included  for  ^  >  6  to 


ensure  asymptotic  exactness  when  C^^are  exact.  The  damp¬ 
ing  functions  fd,  n  for  determining  the  range  of  the  dispersion 
correction  is  given  as 


where  Sr,  n  is  the  adjustable  order-dependent  scaling  param¬ 
eter  of  the  cutoff  radii,  .  In  DFT-D3,  Sr,  s  is  set  to  unity 
while  Sr,  6  is  fitted  for  a  given  functional.  The  an  parameters 
(of  14)  are  global  constants.^^’^^ 

III.  COMPUTATIONAL  DETAILS 

Electronic  structure  calculations  were  done  using  the 
aug-cc-pVTZ^^’^^  basis  set.  For  selected  cases,  calculations 
were  carried  out  using  the  larger  correlation-consistent  va¬ 
lence  quadruple-^  (aug-cc-pVQZ)  basis  sets.  Kohn-Sham^^ 
CAMB3LYP  calculations  were  carried  out  as  previously 


TABLE  1.  HE  and  DFT  based  EFP2  binding  energies  (in  kcal/mol)  for  S22  benchmark  set. 


Systems 

Ref.^ 

EFP2-HF'’ 

EFP2-DFT-D3^ 

DFT-DS*! 

Hydrogen  bonded  complexes 

1,  (NH3)2 

-3.15 

-0.52 

-0.38 

-3.43 

2,  (H20)2 

-5.07 

-5.30 

-4.97 

-5.72 

3,  (Formic  acid)2 

-18.81 

-  12.92 

-13.16 

-21.13 

4,  (Formamide)2 

-16.11 

-  15.43 

-  13.44 

-17.51 

5,  (Uracil)2 

-  20.69 

-  20.23 

-20.51 

-  22.27 

6,  2-Pyridoxine  •  2-aminopyridine 

-  17.00 

-  20.92 

-  15.74 

-18.21 

7,  Adenine  •  thymine 

-  16.74 

-  19.22 

-  17.28 

-  17.75 

MUE 

2.33 

1.88 

1.21 

Predominantly  dispersion  complexes 

8,  (CH4)2 

-0.53 

-  1.48 

-0.84 

-0.48 

9,  (C2H4)2 

-1.48 

-3.45 

-2.02 

-1.73 

10,  C6H6  •  CH4 

-1.45 

-1.87 

-2.49 

-1.59 

11,  (C6H6)2  Stack 

-2.62 

3.24 

2.77 

-2.06 

12,  (Pyrazine)2 

-4.20 

1.35 

-1.75 

-3.65 

13,  (Uracil)2  stack 

-9.74 

-11.98 

-6.74 

-  10.85 

14,  Indole  •  benzene  stack 

-4.59 

7.70 

-1.39 

-3.44 

15,  Adenine  •  thymine  stack 

-11.66 

-12.15 

-11.10 

-11.14 

MUE 

3.72 

2.06 

0.48 

Mixed  complexes 

16,  Ethene  •  ethine  (C2v) 

-1.50 

-1.63 

-1.68 

-1.83 

17,  Benzene  •  H2O  (C^) 

-3.29 

-2.61 

-3.59 

-3.51 

18,  Benzene  •  NH3  (Q) 

-2.32 

0.43 

0.56 

-2.43 

19,  Benzene -HCN  (C,) 

-4.55 

-5.71 

-6.31 

-5.05 

20,  (C6H6)2  T— shape 

-2.71 

-3.58 

-2.51 

-2.82 

21,  Indole  •  benzene  T— shape 

-5.62 

-4.24 

-5.98 

-5.51 

22,  (Phenol)2 

-7.09 

-7.11 

-5.82 

-7.66 

MUE 

1.00 

0.99 

0.28 

MUE  (22) 

2.41 

1.66 

0.65 

^References  14  and  16. 


'’EFP2-HF  generated  with  the  aug-cc-pVTZ  basis  set  using  Edmiston-Ruedenberg  localized  orbital;  MUE  of  2.06  kcal/mol 
(0.9  kcal/mol  after  geometry  optimization)  for  EFP2-HF  using  the  6-31+0*  basis  set  for  electrostatic  and  6-31++G(3df,2p) 
basis  set  for  other  terms.  EFP2-HF  generated  with  the  6-311++G(3df,2p)  basis  set  at  the  MP2/aug-cc-pVTZ  geometries:  —5.1 
kcal/mol  (2,  water  dimer),  —2.4  kcal/mol  (11,  stacked  benzene  dimer),  —3.9  kcal/mol  (17,  benzene-water),  —2.9  kcal/mol  (20, 
T— shape  benzene  dimer). 

‘^Generated  with  CAMB3LYP/aug-cc-pVQZ  using  Boys  localized  orbitals  and  D3  dispersion  correction.  Polarization  terms  for 
formic  acid,  uracil,  2-pyridoxine,  2-aminopyridine  thymine,  indole,  C6H6,  pyrazine,  HCN  were  generated  with  CAMB3LYP/aug- 
cc-pVTZ. 

^CAMB3LYP-D3/aug-cc-pVTZ. 
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Hydrogen  bonded  eomplexes 


1  (NH3)2  2  (H20)2  3  (Formic  Acid)2  4  (Formamide)2 


5  (Uracil)2  6  2-Pyridoxine-2-Aminopyridine  7  Adenine -Thymine 

Predominantly  dispersion  eomplexes 

>-  'I  ^  ^  ^ 


> 

8  (CH,), 


9  lOC^H^  CH^ 


12  (Pyrazine)2  13  (Uraeil)2  staek  14  Indole-Benzene  staek  15  Adenine -Thymine  staek 

Mixed  eomplexes 


—  V 


T- 


16  Ethene-Ethine  17  Benzene-H20  18  Benzene-NHj  19  Benzene-HCN 

^  ‘j  $  ^ 

20  T-shape  21  Indole-Benzene  T-shape  22  (Phenol)2 


EIG.  1 .  Structures  of  S22  data  set. 


described.'^^  Polarizability  tensors  for  all  DFT  based  polariza¬ 
tion  potentials  were  obtained  using  The  finite-field  method, 
except  for  2-Pyridoxine.  Due  to  unphysically  large  polariz¬ 
ability  tensors  obtained  with  the  finite-field  method  for  2- 
Pyridoxine,  its  tensors  were  obtained  using  the  CPKS 
method  with  approximate  response  vectors  that  include  KS 
orbital  energies  and  dipole  moments. 

The  S22^'^’  benchmark  set  reference  data  were  obtained 
by  CCSD(T)/CBS  for  complexes  of  common  molecules 
with  typical  noncovalent  interactions.  The  water-benzene 
clusters  and  WATER27^^  reference  values  were  also  ob¬ 
tained  at  the  CCSD(T)/CBS  level,  except  for  the  clusters  with 
n  =  20,"^^  where  MP2/CBS  was  used.  The  ^-alkane  dimers 
benchmark  values  were  CCSD(T)/CBS  for  ethane  to  butane 
and  a  linear  extrapolation  method  for  longer  ^-alkanes.  The 
structures  of  ethane  to  heptane  were  obtained  at  the  MP2/6- 
31G(d)  level. The  geometries  of  ^-alkane  dimers  were  ob¬ 
tained  at  MP2/6-31  lG(d,p)  level  with  fixed  geometries  of  the 
monomers. All  calculations  were  done  using  a  locally  mod¬ 
ified  version  of  the  GAMES  program. 


IV.  RESULTS  AND  DISCUSSION 
A.  S22  benchmark 

The  S22  benchmark  in  Table  I  includes  results  for  22 
typical  complexes  (Figure  1)  that  can  be  divided  into  three 
subsets,  hydrogen  bonded  complexes,  predominantly  disper¬ 
sion  complexes,  and  mixed  complexes  with  electrostatic  and 
dispersion  interactions.  The  performance  of  different  EFP2 
models  for  binding  energies  in  the  S22  set  are  given  along 
with  the  results  for  the  CAMB3LYP  functional  reported  by 
Goerigk  and  Grimme.^^  The  overall  MUEs  for  22  typical  non¬ 
covalent  complexes  are  given  along  with  the  MUEs  for  three 
subsets.  We  first  discuss  the  overall  performance  of  binding 
energies  and  then  consider  specific  complexes.  An  overall 
MUE  of  2.41  kcal/mol  was  obtained  using  EEP2-HE  with  the 
aug-cc-pVTZ  basis  set.  This  is  comparable  with  the  MUE  of 
2.06  kcal/mol  for  EEP2-HE  using  the  6-3 1+G*  basis  set  for 
electrostatics  and  6-31-h-hG(3df,2p)  basis  set  for  other  terms 
(EEP2-HE-MX).  The  MUE  for  EEP2-HE-MX  was  reduced  to 
0.9  kcal/mol  after  geometry  optimization.  Eor  the  DET  based 
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5  C6H6  (H20)3,  W3Blb 


6  C6H6  (H20)3,  W3B1c  7  (C6H6)2  H20,  WlB2a 


9(C6H6)2(H20)2,W2B2b 


10  (C6H6)2  (H20)2,  W2B2c 


FIG.  2.  Structures  of  water-benzene  clusters. 


EFP2  method,  the  MUE  of  1.66  kcal/mol  is  obtained  for  a 
potential  generated  by  CAMB3LYP/aug-cc-pVQZ  with  Boys 
localized  orbitals  and  the  dispersion  correction  using  the  D3 
(5^,6  =  1.217,  S6  =  1.0,  ^8  =  1.378)  parameters'^  (EFP2- 
CAMB3LYP-D3).  A  MUE  of  0.98  kcal/mol  is  obtained  with 
EFP2-CAMB3LYP-D3  after  geometry  optimization.  The  D3 
attractive  dispersion^^  contributions  to  the  binding  energies 
depend  on  the  type  of  complex  and  are  smaller  than  the  corre¬ 
sponding  HE  values  (see  Tables  IS  and  2S  in  the  supplemen¬ 
tary  material'^'^).  The  errors  for  EFP2-CAMB3LYP-D3  are 
comparable  to  those  of  EFP2-HF-MX.  However,  the  D3  dis¬ 
persion  corrections  were  obtained  at  a  significantly  reduced 
computational  cost  in  comparison  to  the  coupled  perturbed 
HE  method  (CPHE)  for  obtaining  dispersion  coefficients,  ex¬ 
pressed  in  terms  of  the  dynamic  polarizabilties  over  the  imag¬ 
inary  frequency  range,^  which  is  currently  used  for  EEP2-HE. 
The  EFP2  errors  are  higher  than  the  reported^^  MUE  of  0.67 
kcal/mol  for  the  D3  correction  of  CAMB3LYP/def2-QZVP 
(MUE  of  2.52  kcal/mol)  using  the  RI-JK  approximation.  This 
is  quite  small  but  larger  than  MUE  for  the  lowest  BLYP-D3 
(0.23  kcal/mol)^^  functional.  The  MUE  of  0.65  kcal/mol 
is  obtained  for  CAMB3LYP-D3/aug-cc-pVTZ  without  using 
the  RI-JK  approximation. 


B.  Water-benzene  benchmark 

In  this  section,  we  examine  the  water-benzene  clusters 
(Figure  2)  and  related  complexes.  A  number  of  com¬ 
plexes  (water  dimer,  benzene  dimers,  water-benzene)  in 


the  S22  set  have  been  considered  using  EFP2-HF  with 
the  6-311-h-hG(3df,2p)  basis  set  and  MP2/aug-cc-pVTZ 
geometry.  The  reported^^  EFP2-HF  optimized  water 
dimer  classical  binding  energy  of  5.1  kcal/mol  is  in  good 
agreement  with  our  computed  EFP2-HF  (—5.3  kcal/mol) 
and  EFP2-CAMB3LYP-D3  (-4.97  kcal/mol)  values  (Table 
I).  For  the  lowest  T-shape  structure  of  the  benzene  dimer, 
we  obtain  a  binding  energy  (without  zero-point  correc¬ 
tion,  Table  I)  of  -2.51  (-3.6  kcal/mol)  EFP2-CAM-D3 


TABLE  II.  HF  and  DFT  based  EFP2  binding  energies  (in  kcal/mol)  for  the 
water-benzene  clusters  benchmark  set. 


Systems 

Ref.’' 

EFP2-HF*> 

EEP2-HE'^ 

EFP2-DFT-D3'= 

1,  WlBla 

-3.0 

-3.5 

-3.9 

-2.7 

2,  WlBlb 

-  1.1 

-2.8 

-2.6 

-3.2 

3,  W2B1 

-8.9 

-10.7 

-10.9 

-9.5 

4,  W3Bla 

-16.0 

-21.0 

-20.7 

-19.9 

5,  W3Blb 

-17.3 

-20.9 

-19.6 

-19.5 

6,  W3Blc 

-16.7 

-20.4 

-19.4 

-18.9 

7,  WlB2a 

-7.0 

-8.4 

1 

00 

bo 

-8.0 

8,  W2B2a 

-13.0 

-16.7 

-17.3 

-16.0 

9,  W2B2b 

-10.1 

-17.0 

-16.0 

-15.8 

10,  W2B2c 

-11.4 

-13.8 

-14.7 

-12.5 

MUE 

3.1 

2.9 

2.2 

^Reference  17. 

’’Generated  with  the  aug-cc-pVQZ  basis  set  at  the  MP2/cc-pVTZ  geometries. 
‘^Generated  with  the  6-311-1— |-G(3df,2p)  basis  set  at  the  MP2/aug-cc-pVTZ 
geometries.’^ 

^Generated  with  CAMB3LYP/aug-cc-pVQZ  (CAMB3LYP/aug-cc-pVTZ  for  polariza¬ 
tion  potential)  using  Boys  localized  orbitals  and  D3  dispersion  correction. 
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1  (H^O)^ 


2  (H20)3 


5  (H20)6  Prism 


6(H20)6  Cage 


7  (H20)6  Book 


8  (H20)6  Cyclic 


10  (H20)8  Cube  11  (H20)2o  Dodecahedron 


FIG.  3.  Structures  of  water  clusters. 


(EFP2-HF).  EFP2-CAM-D3  predicts  the  stacked  structure 
to  be  unbound  at  the  S22  benchmark  geometry.  How¬ 
ever,  the  EFP2-CAM-D3  (EFP2-HF)  binding  energies  of 
—2.2  kcal/mol  (—1.9  kcal/mol)  and  —2.8  kcal/mol 
(—4.0  kcal/mol)  were  obtained  after  geometry  optimization 
for  the  stacked  and  T- shape  dimers,  respectively.  These  values 
are  in  good  agreement  with  the  previously  reported  EFP2- 
HF  values  of  —2.4  kcal/mol  (stacked)  and  —2.9  kcal/mol 
(T-shape).^^  For  benzene-water  dimer,  the  lowest  EFP2-HF 
classical  binding  energy  of  —3.9  kcal/mol  was  reported 
for  the  geometry  with  water  on  top  of  the  benzene  ring.^^ 
Our  computed  EFP2-HF  values  of  —2.6  kcal/mol  (Ta¬ 
ble  I)  and  —3.5  kcal/mol  (Table  II)  were  obtained  at  the 
MP2/cc-pVTZ  structure  in  the  S22  benchmark  and  at  the 
optimized  EFP2  geometry,  respectively.  The  corresponding 
EFP2-CAMB3FYP-D3  values  are  -3.6  kcal/mol  (Table  I) 
and  —2.7  kcal/mol  (Table  II).  For  complexes  with  multiple 
water  and  benzene  molecules,  their  binding  energies  are 
generally  overestimated  by  both  the  HF  and  DFT  based  EFP2 
methods.  Overall,  MUE  of  about  3  kcal/mol  and  2  kcal/mol 
are  obtained  with  EFP2-HF  and  EFP2-CAMB3FYP-D3 
(Table  II). 


C.  Water  clusters 

In  this  validation  set,  we  considered  the  14  neutral 
(H20)n  (Figure  3)  in  the  so-called  WATER27  benchmark.  The 
benchmark  set  provides  a  rigorous  test  for  hydrogen  bonding 
in  small  (H20)2,  mid- size  (H20)n,  and  the  large  size  of  the 
clusters  up  to  n  =  20.  The  MUEs  for  HF  and  DFT  based 


TABLE  III.  MUEs  for  binding  energies  (in  kcal/mol)  for  the  14  neutral  wa¬ 
ter  clusters  in  the  WATER27  benchmark  set  using  EEP2  and  other  methods.^ 


Basis  set/localization  method 

Edmiston-Ruedenberg 

Boys 

Pipek-Mezey 

EFP2-CAMB3LYP-D3 

aug-cc-pVTZ 

5.8 

9.8 

7.3 

aug-cc-pVQZ 

7.1 

2.0 

15.3 

EFP2-HF 

aug-cc-pVTZ 

10.5 

4.7 

33.7 

aug-cc-pVQZ 

1.8 

11.0 

39.5 

^MUEs  for  other  methods:  11.8  kcal/mol  (EFPl-HF),  43.9  kcal/mol  (EFPl-DFT), 
13.2  kcal/mol  (EFP2-HF  using  the  6-3 U-G*  basis  set  for  electrostatic  with  an  overlap 
based  damping  correction  and  6-31-1— |-G(3df,2p)  basis  set  for  other  terms,  20.1  kcal/mol 
(POL5P),  6.3  kcal/mol,  6.3  kcal/mol  (TIP5P),  5.7  (TIP5PE),  4.2  kcal/mol  (SPC), 
3.3  kcal/mol  (SPCE),  1.9  kcal/mol  (TIP3P),  1.6  kcal/mol  (AMOEBA). 
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TABLE  IV.  EFP2  binding  energies  (in  kcal/mol)  generated  by  HE  and  CAMB3LYP  with  the  aug-cc-pVQZ  basis  set  using  Edmiston-Ruedenberg  and  Boys 
localized  orbitals,  respectively,  for  neutral  water  clusters  in  the  WATER27  benchmark  set  compared  with  other  methods. 


Systems 

Ref.^ 

PW6B95'’ 

KS 

D3 

EFPL 

HE 

DFT 

EFP2 

HE 

EFP2 

CAMB3LYP 

1,  (H20)2 

-5.0 

-4.7 

-5.0 

-5.0 

-7.7 

-4.9 

-5.0 

2,  (H20)3  cyclic 

-15.8 

-14.7 

-15.7 

-13.6 

-22.1 

-15.9 

-15.6 

3,  (H20)4  cyclic 

-27.4 

-25.8 

-27.3 

-24.1 

-39.7 

-26.7 

-25.6 

4,  (H20)5  cyclic 

-35.9 

-34.0 

-36.0 

-31.8 

-52.8 

-35.5 

-36.5 

5,  (H20)6  prism 

-46.0 

-42.0 

-45.9 

-40.9 

-68.7 

-45.9 

-45.1 

6,  (H20)6  cage 

-45.8 

-42.1 

-45.6 

-40.2 

-68.3 

-45.4 

-45.6 

7,  (H20)6  book 

-45.3 

-42.2 

-45.2 

-40.2 

-67.3 

-44.5 

-43.5 

8,  (H20)6  cyclic 

-44.3 

-42.2 

-44.6 

-39.7 

-65.2 

-43.5 

-45.7 

9,  (H20)8  cube  {P2d) 

-72.6 

-66.0 

-71.6 

-64.1 

-  109.3 

-71.3 

-69.9 

10,  (H20)8  cube  {S4) 

-72.6 

-65.9 

-71.6 

-64.1 

-  109.2 

-71.3 

-70.5 

11,  (H20)2o  dodecahedron 

-200.1 

-182.4 

-  197.6 

-  176.8 

-301.6 

-  197.5 

-  203.5 

12,  (H20)2o  fused  cubes 

-212.6 

-184.5 

-  203.6 

-  182.6 

-318.5 

-  208.3 

-  207.7 

13,  (H20)2o  face— sharing 

-215.0 

-185.5 

-  204.3 

-  183.6 

-  320.2 

-211.2 

-212.4 

14,  (H20)2o  edge— sharing 

-217.9 

-187.8 

-  206.2 

-  184.8 

-  320.7 

-  208.8 

-211.9 

MUE 

9.8 

2.6 

11.8 

43.9 

1.8 

2.0 

‘^Reference  15. 

’’Reference  13. 

‘^EFPl-HF  results  using  HF/DH(d,p)//HF/DH(d,p):  -42.4  (5,  prism),  -41.9  (6,  cage),  -41.5  (7,  book),  -41.1  (8,  cyclic),  -66.3  (9,  D2d  cube),  -66.3  (10,  S4  cube).^^ 


EFP2  potentials  using  different  localized  orbitals  and  basis 
sets  (aug-cc-pVXZ,  X  =  T,  Q)  are  listed  in  Table  III,  along 
with  the  errors  obtained  with  other  methods.  For  a  given  ba¬ 
sis  set,  the  variations  of  the  results  among  different  localiza¬ 
tion  schemes  can  be  attributed  to  the  polarization  potential 
(Eq.  (10)),  that  is  expanded  in  terms  of  EMO  polarizability 
tensors.  Although  polarization  energies  obtained  with  EMO 
polarizabilities  were  found  to  be  more  accurate  than  the  en¬ 
ergies  expanded  with  polarizabilties  at  the  center  of  mass,^ 
different  LMOs  result  in  different  polarizability  tensors  and 
resulting  energies.  Similarly,  the  dispersion  potential  obtained 
with  EMO  dynamic  polarizabilities,^  is  an  additional  con¬ 
tributing  factor  of  the  deviations  for  the  HE  based  EFP2.  Since 
both  static  and  dynamic  polarizabilities  are  basis  set  sensitive, 
the  resulting  energies  are  basis  dependent. 

The  errors  for  water  clusters  are  generally  improved 
upon  going  from  the  aug-cc-pVTZ  to  aug-cc-pVQZ  basis 
set,  except  for  values  obtained  with  the  Pipek-Mezey  local¬ 
ized  orbitals  and  the  Boys  localized  orbitals  for  EFP2-HF. 
The  lowest  MUEs  of  1.8  kcal/mol  and  2.0  kcal/mol  were  ob¬ 
tained  with  the  Edminston-Ruedenberg  and  Boys  localized 
aug-cc-pVQZ  orbitals  for  HE  and  CAMB3LYP  based  EFP2 
potentials,  respectively.  These  errors  compare  well  to  the 
best  estimates  from  DFT-D3  (2.6  kcal/mol  for  PW6B95,  see 
Table  IV)  and  other  (1.6  kcal/mol  for  AMOEBA)  potentials. 
The  binding  energies  for  EFP2  obtained  with  HE  and  the 
CAMB3LYP  functional  are  listed  in  Table  IV  along  with  the 
results  for  EFPl,  PW6B95  functional,  and  reference  values. 
We  obtain  a  value  of  11.8  kcal/mol  for  EFPl -HE,  which  has 
been  applied  to  study  (H20)n  for  n  =  6-20."^^  The  MUE 
(10.5  kcal/mol,  see  Table  III)  for  a  EFP2-HF  potential  using 
the  aug-cc-pVTZ  basis  set  and  Edminston-Ruedenberg  local¬ 
ized  orbitals  is  slightly  reduced  and  is  comparable  to  the  error 
(9.8  kcal/mol,  see  Table  IV)  for  PW6B95  without  dispersion 
corrections.  With  the  D3  dispersion  correction,  the  MUE  for 


PW6B95  reduces  to  2.6  kcal/mol.  EFP2-CAM-D3  delivers 
outstanding  improvements  over  EFP2-CPHF  and  EFPl  fitted 
potentials  obtained  for  the  water  clusters  benchmark  set.  The 
MUE  of  2.0  kcal/mol  for  EFP2-CAMB3LYP  is  smaller  than 
the  corresponding  value  (2.6  kcal/mol)  obtained  for  the  most 
accurate  PW6B95-D3  functional  in  the  WATER27  bench¬ 
mark  set.  These  errors  are  quite  small,  considering  the  average 
binding  energy  of  —89.7  kcal/mol.  The  poor  EFP1-B3LYP  re¬ 
sults  might  be  due  to  known  deficiencies  of  DFT  densities  at 
large  interatomic  distances  in  complexes  for  functionals  with 
asymptotically  incorrect  exchange-correlation  potentials. 

The  ability  to  predict  relative  energies  of  low-energy  iso¬ 
mers  is  important.  Thus,  we  examine  the  performance  of  dif¬ 
ferent  models  for  relative  stability  of  (H20)6  and  (H20)2o 
conformers.  For  (H20)6,  CCSD(T)  gives  the  following  energy 
ordering  for  the  four  isomers:  prism  <  cage  <  book  <  cyclic 


1  (Ethane)2 


2  (Propane)2 


3  (Butane)2 


5  (Hexane)2 


6  (Heptane)2 


FIG.  4.  Structures  of  ^-alkane  dimers. 
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TABLE  V.  HE  and  DFT  based  EFP2  binding  energies  (in  kcal/mol)  using 
the  aug-cc-pVTZ  basis  set  for  the  n-alkane  dimers.  Values  in  parentheses  are 
obtained  using  the  aug-cc-pVQZ  basis  set. 


Systems 

Ref."* 

EFP2-HF^ 

EFP2-DFT-D3'= 

1,  Ethane 

-1.22 

-0.16  (-0.28) 

-  0.55  (-0.70) 

2,  Propane 

-  1.87 

-2.51  (-2.43) 

-  2.27  (-2.20) 

3,  Butane 

-2.74 

-5.21 

-3.33  (-3.41) 

4,  Pentane 

-3.80 

-6.28 

-4.21 

5,  Hexane 

-4.66 

-8.66 

-5.32 

6,  Heptane 

-5.52 

-9.95 

-5.66 

MUE 

2.51 

0.48 

^Reference  18. 

'’Generated  with  the  aug-cc-pVTZ  basis  set  at  the  MP2/cc-pVTZ  geometries  using 
Edmiston-Ruedenberg  localized  orbitals. 

‘^Generated  with  CAMB3LYP/aug-cc-pVTZ  using  Boys  localized  orbitals  and  D3  dis¬ 
persion  correction. 

(5-8  in  Figure  3).  The  EFPl  and  EFP2-HF  models  correctly 
predict  the  ordering.  However,  EFPl -DFT  excessively  overes¬ 
timates  the  binding  energies  while  EFPl -HE  underestimates 
the  energies.  EFP2-CAMB3LYP-D3  predicts  small  energy 
gaps  among  the  isomers.  For  (H20)2o,  the  reference  bind¬ 
ing  energies  predicted  by  MP2/CBS"^^  are  available  for  the 
stable  isomers  of  the  four  different  topologies.  EFP2  models 
predict  the  face- sharing  instead  of  the  edge- sharing  pentago¬ 
nal  prisms  as  the  lowest  energy  isomer.  EFP2-CAMB3LYP- 
D3  and  EFP2-HF  models  underestimate  the  binding  energies 
of  the  lowest  energy  (H20)2o  isomer  by  6.0  kcal/mol  and 
9.1  kcal/mol,  respectively. 

D.  n-Alkane  dimers  benchmark 

The  ^-alkane  dimers  benchmark^  ^  set  contains  six 
dimerization  energies  for  all-fra/t5’-/2-alkane  molecules  from 
ethane  to  heptane  (Figure  4).  The  CCSD(T)/CBS  reference 
energies  were  MP2/CBS  values  corrected  by  the  (CCSD(T) 
—  MP2)/cc-pVTZ  estimate  for  ethane  to  butane  and  a  linear 
extrapolation  of  the  coupled  cluster  correction  (kcal/mol)  by 
—0.068  -h  0.0616  V,  where  x  is  the  number  of  carbon  atoms 
in  the  monomer,  for  the  longer  chains.  These  reference  val¬ 
ues  are  listed  in  Table  V  along  with  the  EFP2  results.  Inter¬ 
action  energies  for  these  complexes  vary  from  —1.3  kcal/mol 
to  —5.6  kcal/mol.  The  EFP2-CAMB3LYP-D3  interaction  en¬ 
ergies  are  in  good  agreement  with  the  reference  values  in 
Table  V.  We  found  that  EFP2  interaction  energies  do  not 
vary  significantly  with  basis  sets  and  localization  (Boys  and 
Emiston-Ruedenberg)  methods.  The  EEP2/CAMB3LYP  po¬ 
tential  generated  with  the  aug-cc-pVTZ  basis  set  has  the  MUE 
of  0.5  kcal/mol.  This  is  much  smaller  than  the  MUE  of  2.5 
kcal/mol  obtained  with  the  EEP2/HE  potential  using  the  same 
basis  set. 

V.  SUMMARY  AND  CONCLUSIONS 

A  method  for  generating  effective  fragment  potentials 
based  on  density  functional  theory  (EEP2-DET)  is  developed 
for  the  treatment  of  solvent  effects.  EEP2-DET  includes  elec¬ 
trostatic,  exchange-repulsion,  polarization,  and  dispersion. 


which  are  generated  for  a  chosen  DET  functional  for  a  given 
isolated  molecule  with  a  computational  cost  comparable  to 
that  of  the  HE  based  EEP2  method.  In  analogy  to  the  HE  based 
EEP2,  the  electrostatic  potential  is  modeled  with  multipolar 
expansion  at  each  atomic  center  and  bond  midpoint  by  Stone's 
DET  distributed  multipolar  analysis.  The  exchange-repulsion 
potential  between  two  fragments  is  composed  of  the  overlap 
and  kinetic  energy  integrals  and  the  nondiagonal  KS  matri¬ 
ces  in  the  localized  molecular  orbital  basis.  The  polarization 
potential  is  derived  from  the  static  molecular  polarizability 
from  DET  finite  field  calculations  and  approximate  coupled 
perturbed  Kohn-Sham  equations.  The  dispersion  potential  in¬ 
corporates  the  low  cost  D3  dispersion  correction  of  Grimme 
et  that  extends  the  range  applicability  of  EEP  to  larger 
system  sizes  and  provides  comparable  or  improved  results  for 
EEP2-DET  over  previous  HP  based  EPP2  method,  which  has 
MUEs  of  2.41,  3.1,  1.8,  and  2.5  kcal/mol  for  the  S22,  water- 
benzene  clusters,  water  clusters,  and  ^-alkane  dimers  bench¬ 
mark  sets,  respectively.  A  DPT  based  EPP2  potential  using  the 
CAM-B3LYP  functional  was  evaluated  with  the  same  bench¬ 
mark  sets  and  yielded  the  corresponding  MUEs  of  1.7,  2.2, 
2.0,  and  0.5  kcal/mol  for  the  respective  sets. 
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